rm(list = ls())
source("Accepted Paper Replication Archive/code/utils.R")

## Load and clean Qualtrics survey data (main sample)
qualtrics <- read_csv("Accepted Paper Replication Archive/clean_data/qualtrics.csv") %>%
  as.data.frame()
dim(qualtrics)

table(qualtrics$trap, exclude = NULL) ## check whether speed trap worked
length(qualtrics$trap[qualtrics$trap == "Extremely interested,Very interested"]) == nrow(qualtrics)

###
### code respondent attributes
###

table(qualtrics$age, exclude = NULL)
qualtrics$age2 <- as.numeric(as.character(qualtrics$age))
qualtrics$age2[qualtrics$age == "100 or over"] <- 100
summary(qualtrics$age2)

age <- median(qualtrics$age2)

table(qualtrics$gender, exclude = NULL)
qualtrics$female <- NA
qualtrics$female[qualtrics$gender == "Female"] <- 1
qualtrics$female[qualtrics$gender == "Male"] <- 0
table(qualtrics$female, exclude = NULL)

female <- mean(qualtrics$female)
female

table(qualtrics$race, exclude = NULL)

qualtrics$hispanic <- NA
qualtrics$hispanic[!is.na(qualtrics$race)] <- 0
qualtrics$hispanic[grep("Hispanic", qualtrics$race)] <- 1
table(qualtrics$hispanic, exclude = NULL)
hispanic <- mean(qualtrics$hispanic)

qualtrics$nhwhite <- NA
qualtrics$nhwhite[!is.na(qualtrics$race)] <- 0
qualtrics$nhwhite[qualtrics$race == "White/Caucasian"] <- 1
qualtrics$nhwhite[qualtrics$hispanic == 1] <- 0
table(qualtrics$nhwhite, exclude = NULL)
nhwhite <- mean(qualtrics$nhwhite)

qualtrics$nhblack <- NA
qualtrics$nhblack[!is.na(qualtrics$race)] <- 0
qualtrics$nhblack[qualtrics$race == "African American"] <- 1
qualtrics$nhblack[qualtrics$hispanic == 1] <- 0
table(qualtrics$nhblack, exclude = NULL)
nhblack <- mean(qualtrics$nhblack)

qualtrics$nhasian <- NA
qualtrics$nhasian[!is.na(qualtrics$race)] <- 0
qualtrics$nhasian[qualtrics$race == "Asian"] <- 1
qualtrics$nhasian[qualtrics$hispanic == 1] <- 0
table(qualtrics$nhasian, exclude = NULL)
nhasian <- mean(qualtrics$nhasian)

qualtrics$other <- NA
qualtrics$other <- as.numeric(I(qualtrics$hispanic == 0 & qualtrics$nhwhite == 0 & qualtrics$nhblack == 0 & qualtrics$nhasian == 0))
table(qualtrics$other, exclude = NULL)
other <- mean(qualtrics$other)


# check for mutual exclusivity
sum(
  table(qualtrics$hispanic)[2], table(qualtrics$nhwhite)[2],
  table(qualtrics$nhblack)[2], table(qualtrics$nhasian)[2],
  table(qualtrics$other)[2]
) == nrow(qualtrics[!is.na(qualtrics$race), ])

table(qualtrics$education, exclude = NULL)
qualtrics$hs <- 1
qualtrics$hs[qualtrics$education == "Less than High School"] <- 0
table(qualtrics$hs)

table(qualtrics$education, exclude=NULL)#no missings
qualtrics$ba <- 1
qualtrics$ba[qualtrics$education == "High School / GED"] <- 0
qualtrics$ba[qualtrics$education == "Less than High School"] <- 0
qualtrics$ba[qualtrics$education == "Some College"] <- 0
qualtrics$ba[qualtrics$education == "2-year College Degree"] <- 0
table(qualtrics$ba)
ba <- mean(qualtrics$ba)

table(qualtrics$ideo, exclude = NULL)
qualtrics$ideo2 <- NA
qualtrics$ideo2[qualtrics$ideo == "Very conservative"] <- 1
qualtrics$ideo2[qualtrics$ideo == "Conservative"] <- 2
qualtrics$ideo2[qualtrics$ideo == "Moderate"] <- 3
qualtrics$ideo2[qualtrics$ideo == "Liberal"] <- 4
qualtrics$ideo2[qualtrics$ideo == "Very liberal"] <- 5
table(qualtrics$ideo2, exclude = NULL)

## dichotomous party measure
table(qualtrics$party, exclude = NULL)
table(qualtrics$closer, exclude = NULL)
qualtrics$party2 <- NA
qualtrics$party2[qualtrics$party == "Democrat"] <- "dem"
qualtrics$party2[qualtrics$party == "Republican"] <- "rep"
qualtrics$party2[qualtrics$closer == "Democratic Party"] <- "dem"
qualtrics$party2[qualtrics$closer == "Republican Party"] <- "rep"
qualtrics$party2[qualtrics$closer == "Neither"] <- "ind"
table(qualtrics$party2, exclude = NULL)

dem <- mean(I(qualtrics$party == "Democrat"))
rep <- mean(I(qualtrics$party == "Republican"))
ind <- mean(I(qualtrics$party == "Independent"))

dem2 <- mean(I(qualtrics$party2 == "dem"))
rep2 <- mean(I(qualtrics$party2 == "rep"))
ind2 <- mean(I(qualtrics$party2 == "ind"))

## Table B1 ----------------------------------------


vars <-
  list(
    female,
    age,
    ba,
    hispanic,
    nhwhite,
    nhblack,
    nhasian,
    other,
    dem2,
    rep2,
    ind2
  )
demo <- as.data.frame(matrix(nrow = length(vars), ncol = 3))
colnames(demo) <- c("Variable", "Sample", "U.S.")
demo$Variable <-
  c(
    "Female",
    "Age (median years)",
    "At Least BA",
    "Hispanic",
    "Non-Hispanic White",
    "Non-Hispanic Black",
    "Non-Hispanic Asian",
    "Other Race",
    "Democrat",
    "Republican",
    "Independent"
  )
for (i in 1:length(vars)) {
  if (i == 2) {
    demo[i, "Sample"] <- paste0(round(vars[[i]], digits = 2), "")
  }

  if (i != 2) {
    demo[i, "Sample"] <- paste0(round(vars[[i]], digits = 2) * 100, "%")
  }
}

demo

## party ID in 2020 ANES: https:/electionstudies.org/resources/anes-guide/top-tables/?id=21

anes <- read_csv("Accepted Paper Replication Archive/clean_data/anes.csv") %>%
  as.data.frame()

# 7-point party ID
anes$party7 <- as.numeric(anes$V201231x)
table(anes$party7)

mean_dems_anes <- mean(I(anes$party7 %in% c(1, 2, 3)))
mean_reps_anes <- mean(I(anes$party7 %in% c(5, 6, 7)))
mean_ind_anes <- mean(I(anes$party7 %in% c(4)))

demo[demo$Variable == "Democrat", "U.S."] <-
  paste0(round(mean_dems_anes, digits = 2) * 100, "%")
demo[demo$Variable == "Republican", "U.S."] <-
  paste0(round(mean_reps_anes, digits = 2) * 100, "%")
demo[demo$Variable == "Independent", "U.S."] <-
  paste0(round(mean_ind_anes, digits = 2) * 100, "%")

# 2019 1-year ACS estimates
census_2019_1year_acs <- read_csv("Accepted Paper Replication Archive/clean_data/census_2019_1year_acs.csv") %>%
  as.data.frame()

nhwhite_count <-
  as.numeric(census_2019_1year_acs[census_2019_1year_acs$POPGROUP_LABEL == "White alone, not Hispanic or Latino", "S0201_001E"])
nhblack_count <-
  as.numeric(census_2019_1year_acs[census_2019_1year_acs$POPGROUP_LABEL == "Black or African American alone, not Hispanic or Latino", "S0201_001E"])
nhasian_count <-
  as.numeric(census_2019_1year_acs[census_2019_1year_acs$POPGROUP_LABEL == "Asian alone, not Hispanic or Latino", "S0201_001E"])
hispanic_count <-
  as.numeric(census_2019_1year_acs[census_2019_1year_acs$POPGROUP_LABEL == "Hispanic or Latino (of any race) (200-299)", "S0201_001E"])
total <-
  as.numeric(census_2019_1year_acs[census_2019_1year_acs$POPGROUP_LABEL == "Total population", "S0201_001E"])

demo[demo$Variable == "Hispanic", "U.S."] <-
  paste0(round(hispanic_count / total, digits = 2) * 100, "%")
demo[demo$Variable == "Non-Hispanic White", "U.S."] <-
  paste0(round(nhwhite_count / total, digits = 2) * 100, "%")
demo[demo$Variable == "Non-Hispanic Black", "U.S."] <-
  paste0(round(nhblack_count / total, digits = 2) * 100, "%")
demo[demo$Variable == "Non-Hispanic Asian", "U.S."] <-
  paste0(round(nhasian_count / total, digits = 2) * 100, "%")
demo[demo$Variable == "Other Race", "U.S."] <-
  paste0(
    round(
      1 - hispanic_count / total - nhwhite_count / total - nhblack_count / total - nhasian_count /
        total,
      digits = 2
    ) * 100,
    "%"
  )



## other census_2019_1year_acssus variables: https://www.census_2019_1year_acssus.gov/quickfacts/fact/table/US/PST045221
# has % with BA and gender
census_ba_gender <- read_csv("Accepted Paper Replication Archive/clean_data/census_ba_gender.csv") %>%
  as.data.frame()
female.us <- census_ba_gender[
  census_ba_gender$Fact == "Female persons, percent",
  "United.States"
][1]
female.us <- gsub("%", "", female.us)
female.us <- round(as.numeric(female.us), 0)

ba.us <- census_ba_gender[
  census_ba_gender$Fact == "Bachelor's degree or higher, percent of persons age 25 years+, 2018-2022",
  "United.States"
][1]
ba.us <- gsub("%", "", ba.us)
ba.us <- round(as.numeric(ba.us), 0)

demo[demo$Variable == "At Least BA", "U.S."] <- paste0(ba.us, "%")
demo[demo$Variable == "Female", "U.S."] <- paste0(female.us, "%")

# get median age
census_median_age <- read_csv("Accepted Paper Replication Archive/clean_data/census_median_age.csv") %>%
  as.data.frame()

med_age <- as.numeric(census_median_age$S0101_C01_032E[2])
demo[demo$Variable == "Age (median years)", "U.S."] <- round(med_age, digits = 0)
demo

print(xtable(demo), include.rownames = FALSE)

###
## code TREATMENTS in Dose Response Experiment
###

# dose for DOD
table(qualtrics$pct_dose_dod, exclude = NULL)
qualtrics$pct_dose_dod[is.na(qualtrics$pct_dose_dod)] <- "control"
qualtrics$pct_dose_dod <-
  factor(qualtrics$pct_dose_dod,
    levels = c("control", "20", "30", "40", "50", "60", "70", "80")
  )
table(qualtrics$pct_dose_dod)

# dose for Treasury
table(qualtrics$pct_dose_treas, exclude = NULL)
qualtrics$pct_dose_treas <- as.character(qualtrics$pct_dose_treas)
qualtrics$pct_dose_treas[is.na(qualtrics$pct_dose_treas)] <- "control"
qualtrics$pct_dose_treas <-
  factor(qualtrics$pct_dose_treas,
    levels = c("control", "20", "30", "40", "50", "60", "70", "80")
  )
table(qualtrics$pct_dose_treas)

# dose for HHS
table(qualtrics$pct_dose_hhs, exclude = NULL)
qualtrics$pct_dose_hhs <- as.character(qualtrics$pct_dose_hhs)

qualtrics$pct_dose_hhs[is.na(qualtrics$pct_dose_hhs)] <- "control"
qualtrics$pct_dose_hhs <-
  factor(qualtrics$pct_dose_hhs,
    levels = c("control", "20", "30", "40", "50", "60", "70", "80")
  )
table(qualtrics$pct_dose_hhs)

# dose for DOE
table(qualtrics$pct_dose_educ, exclude = NULL)
qualtrics$pct_dose_educ <- as.character(qualtrics$pct_dose_educ)
qualtrics$pct_dose_educ[is.na(qualtrics$pct_dose_educ)] <- "control"
qualtrics$pct_dose_educ <-
  factor(qualtrics$pct_dose_educ,
    levels = c("control", "20", "30", "40", "50", "60", "70", "80")
  )
table(qualtrics$pct_dose_educ)

# gender treatment for DOD
table(qualtrics$exp_1_gender_dod) ## gender % emphasized
qualtrics$exp_1_gender_dod2 <- NA
qualtrics$exp_1_gender_dod2[qualtrics$exp_1_gender_dod == "men"] <- "men"
qualtrics$exp_1_gender_dod2[qualtrics$exp_1_gender_dod == "women"] <- "women"
qualtrics$exp_1_gender_dod2 <-
  factor(qualtrics$exp_1_gender_dod2, levels = c("men", "women"))
class(qualtrics$exp_1_gender_dod2)
table(qualtrics$exp_1_gender_dod2, exclude = NULL)

# gender treatment for Treasury
table(qualtrics$exp_1_gender_treas) ## gender % emphasized
qualtrics$exp_1_gender_treas2 <- NA
qualtrics$exp_1_gender_treas2[qualtrics$exp_1_gender_treas == "men"] <- "men"
qualtrics$exp_1_gender_treas2[qualtrics$exp_1_gender_treas == "women"] <- "women"
qualtrics$exp_1_gender_treas2 <-
  factor(qualtrics$exp_1_gender_treas2, levels = c("men", "women"))
class(qualtrics$exp_1_gender_treas2)
table(qualtrics$exp_1_gender_treas2, exclude = NULL)

# gender treatment for HHS
table(qualtrics$exp_1_gender_hhs) ## gender % emphasized
qualtrics$exp_1_gender_hhs2 <- NA
qualtrics$exp_1_gender_hhs2[qualtrics$exp_1_gender_hhs == "men"] <- "men"
qualtrics$exp_1_gender_hhs2[qualtrics$exp_1_gender_hhs == "women"] <- "women"
qualtrics$exp_1_gender_hhs2 <-
  factor(qualtrics$exp_1_gender_hhs2, levels = c("men", "women"))
class(qualtrics$exp_1_gender_hhs2)
table(qualtrics$exp_1_gender_hhs2, exclude = NULL)

# gender treatment for DOE
table(qualtrics$exp_1_gender_educ) ## gender % emphasized
qualtrics$exp_1_gender_educ2 <- NA
qualtrics$exp_1_gender_educ2[qualtrics$exp_1_gender_educ == "men"] <- "men"
qualtrics$exp_1_gender_educ2[qualtrics$exp_1_gender_educ == "women"] <- "women"
qualtrics$exp_1_gender_educ2 <-
  factor(qualtrics$exp_1_gender_educ2, levels = c("men", "women"))
class(qualtrics$exp_1_gender_educ2)
table(qualtrics$exp_1_gender_educ2, exclude = NULL)



## order of 4 agency questions
table(qualtrics$FL_240_DO)
qualtrics$ag_order1 <- NA
qualtrics$ag_order1[substr(qualtrics$FL_240_DO, start = 1, stop = 11) == "exp_1_educ|"] <-
  "educ"
qualtrics$ag_order1[substr(qualtrics$FL_240_DO, start = 1, stop = 10) == "exp_1_hhs|"] <-
  "hhs"
qualtrics$ag_order1[substr(qualtrics$FL_240_DO, start = 1, stop = 10) == "exp_1_dod|"] <-
  "dod"
qualtrics$ag_order1[substr(qualtrics$FL_240_DO, start = 1, stop = 15) == "exp_1_treasury|"] <-
  "treasury"
table(qualtrics$ag_order1, exclude = NULL)

qualtrics$ag_order2 <- NA
qualtrics$ag_order3 <- NA
qualtrics$ag_order4 <- NA

temp <- strsplit(qualtrics$FL_240_DO, "\\|")
head(temp)
for (i in 2:4) {
  for (j in 1:length(temp)) {
    if (temp[[j]][i] == "exp_1_treasury") {
      qualtrics[j, paste0("ag_order", i)] <- "treasury"
    }
    if (temp[[j]][i] == "exp_1_dod") {
      qualtrics[j, paste0("ag_order", i)] <- "dod"
    }
    if (temp[[j]][i] == "exp_1_hhs") {
      qualtrics[j, paste0("ag_order", i)] <- "hhs"
    }
    if (temp[[j]][i] == "exp_1_educ") {
      qualtrics[j, paste0("ag_order", i)] <- "educ"
    }
  }
}

table(qualtrics$ag_order1, exclude = NULL) # agency displayed first
table(qualtrics$ag_order2, exclude = NULL) # agency displayed second
table(qualtrics$ag_order3, exclude = NULL) # agency displayed third
table(qualtrics$ag_order4, exclude = NULL) # agency displayed fourth

# rank treatment
table(qualtrics$rank)
qualtrics$rank <- factor(qualtrics$rank, levels = c("jobs", "top jobs"))
table(qualtrics$rank)

## Outcomes for Dose Response Experiment

# DOD run agency well
table(qualtrics$dod_ability_dept, exclude = NULL)
qualtrics$dod_ability_dept2 <- NA
qualtrics$dod_ability_dept2[qualtrics$dod_ability_dept == "No confidence at all"] <- 1
qualtrics$dod_ability_dept2[qualtrics$dod_ability_dept == "Not too much confidence"] <-
  2
qualtrics$dod_ability_dept2[qualtrics$dod_ability_dept == "Some confidence"] <- 3
qualtrics$dod_ability_dept2[qualtrics$dod_ability_dept == "A lot of confidence"] <- 4
table(qualtrics$dod_ability_dept2, exclude = NULL)

# DOD represent women
table(qualtrics$dod_rep_women, exclude = NULL)
qualtrics$dod_rep_women2 <- NA
qualtrics$dod_rep_women2[qualtrics$dod_rep_women == "No confidence at all"] <- 1
qualtrics$dod_rep_women2[qualtrics$dod_rep_women == "Not too much confidence"] <- 2
qualtrics$dod_rep_women2[qualtrics$dod_rep_women == "Some confidence"] <- 3
qualtrics$dod_rep_women2[qualtrics$dod_rep_women == "A lot of confidence"] <- 4
table(qualtrics$dod_rep_women2, exclude = NULL)

# Treasury run agency well
table(qualtrics$treasury_abil_dep, exclude = NULL)
qualtrics$treasury_abil_dep2 <- NA
qualtrics$treasury_abil_dep2[qualtrics$treasury_abil_dep == "No confidence at all"] <-
  1
qualtrics$treasury_abil_dep2[qualtrics$treasury_abil_dep == "Not too much confidence"] <-
  2
qualtrics$treasury_abil_dep2[qualtrics$treasury_abil_dep == "Some confidence"] <- 3
qualtrics$treasury_abil_dep2[qualtrics$treasury_abil_dep == "A lot of confidence"] <-
  4
table(qualtrics$treasury_abil_dep2, exclude = NULL)

# Treasury represent women
table(qualtrics$treasury_rep_women, exclude = NULL)
qualtrics$treasury_rep_women2 <- NA
qualtrics$treasury_rep_women2[qualtrics$treasury_rep_women == "No confidence at all"] <-
  1
qualtrics$treasury_rep_women2[qualtrics$treasury_rep_women == "Not too much confidence"] <-
  2
qualtrics$treasury_rep_women2[qualtrics$treasury_rep_women == "Some confidence"] <- 3
qualtrics$treasury_rep_women2[qualtrics$treasury_rep_women == "A lot of confidence"] <-
  4
table(qualtrics$treasury_rep_women2, exclude = NULL)

# HHS run agency well
table(qualtrics$hhs_ability_dept, exclude = NULL)
qualtrics$hhs_ability_dept2 <- NA
qualtrics$hhs_ability_dept2[qualtrics$hhs_ability_dept == "No confidence at all"] <- 1
qualtrics$hhs_ability_dept2[qualtrics$hhs_ability_dept == "Not too much confidence"] <-
  2
qualtrics$hhs_ability_dept2[qualtrics$hhs_ability_dept == "Some confidence"] <- 3
qualtrics$hhs_ability_dept2[qualtrics$hhs_ability_dept == "A lot of confidence"] <- 4
table(qualtrics$hhs_ability_dept2, exclude = NULL)

# HHS represent women
table(qualtrics$hhs_rep_women, exclude = NULL)
qualtrics$hhs_rep_women2 <- NA
qualtrics$hhs_rep_women2[qualtrics$hhs_rep_women == "No confidence at all"] <- 1
qualtrics$hhs_rep_women2[qualtrics$hhs_rep_women == "Not too much confidence"] <- 2
qualtrics$hhs_rep_women2[qualtrics$hhs_rep_women == "Some confidence"] <- 3
qualtrics$hhs_rep_women2[qualtrics$hhs_rep_women == "A lot of confidence"] <- 4
table(qualtrics$hhs_rep_women2, exclude = NULL)

# DOE run agency well
table(qualtrics$educ_ability_dept, exclude = NULL)
qualtrics$educ_ability_dept2 <- NA
qualtrics$educ_ability_dept2[qualtrics$educ_ability_dept == "No confidence at all"] <-
  1
qualtrics$educ_ability_dept2[qualtrics$educ_ability_dept == "Not too much confidence"] <-
  2
qualtrics$educ_ability_dept2[qualtrics$educ_ability_dept == "Some confidence"] <- 3
qualtrics$educ_ability_dept2[qualtrics$educ_ability_dept == "A lot of confidence"] <-
  4
table(qualtrics$educ_ability_dept2, exclude = NULL)

# DOE represent women
table(qualtrics$educ_rep_women, exclude = NULL)
qualtrics$educ_rep_women2 <- NA
qualtrics$educ_rep_women2[qualtrics$educ_rep_women == "No confidence at all"] <- 1
qualtrics$educ_rep_women2[qualtrics$educ_rep_women == "Not too much confidence"] <- 2
qualtrics$educ_rep_women2[qualtrics$educ_rep_women == "Some confidence"] <- 3
qualtrics$educ_rep_women2[qualtrics$educ_rep_women == "A lot of confidence"] <- 4
table(qualtrics$educ_rep_women2, exclude = NULL)


### order of experiments

## FL_133 = dose response
## FL_178 = job candidate (nominee)
## FL_169 = charts

# indicator for which experiment was presented first.
qualtrics$order_exp <- NA
qualtrics$order_exp[qualtrics$FL_134_DO == "FL_133|FL_169|FL_178" |
  qualtrics$FL_134_DO == "FL_133|FL_178|FL_169"] <- "dose"
qualtrics$order_exp[qualtrics$FL_134_DO == "FL_178|FL_133|FL_169" |
  qualtrics$FL_134_DO == "FL_178|FL_169|FL_133"] <- "job"
qualtrics$order_exp[qualtrics$FL_134_DO == "FL_169|FL_133|FL_178" |
  qualtrics$FL_134_DO == "FL_169|FL_178|FL_133"] <- "charts"
table(qualtrics$order_exp)

qualtrics$order_exp <- factor(qualtrics$order_exp, levels = c("dose", "job", "charts"))
table(qualtrics$order_exp)

### recode DVs to range from zero to 1 for analysis
qualtrics$dod_ability_dept2 <- recode_0_1(qualtrics$dod_ability_dept2)
qualtrics$dod_rep_women2 <- recode_0_1(qualtrics$dod_rep_women2)
qualtrics$treasury_abil_dep2 <- recode_0_1(qualtrics$treasury_abil_dep2)
qualtrics$treasury_rep_women2 <- recode_0_1(qualtrics$treasury_rep_women2)
qualtrics$hhs_ability_dept2 <- recode_0_1(qualtrics$hhs_ability_dept2)
qualtrics$hhs_rep_women2 <- recode_0_1(qualtrics$hhs_rep_women2)
qualtrics$educ_ability_dept2 <- recode_0_1(qualtrics$educ_ability_dept2)
qualtrics$educ_rep_women2 <- recode_0_1(qualtrics$educ_rep_women2)


## Reshape to long data for analysis of Dose Response Experiment

qualtrics$ID <- 1:nrow(qualtrics) ## respondent ID

temp0 <- matrix(nrow = 4, ncol = 6)
temp0 <- as.data.frame(temp0)
colnames(temp0) <-
  c(
    "ID",
    "agency",
    "gender_treatment",
    "dose",
    "ability",
    "rep_women"
  )
temp0$agency <- c("dod", "treasury", "hhs", "educ")
temp0$ID <- 0
temp0$ag_order <- NA
IDs <- as.numeric(unique(qualtrics$ID))

for (i in 1:length(IDs)) {
  temp <- as.data.frame(matrix(nrow = 4, ncol = 6))
  head(temp)
  colnames(temp) <-
    c(
      "ID",
      "agency",
      "gender_treatment",
      "dose",
      "ability",
      "rep_women"
    )
  temp$agency <- c("dod", "treasury", "hhs", "educ")
  temp$ID <- rep(IDs[i], 4) # 4 observations per respondent
  temp$ag_order <- NA

  # treatments
  temp$gender_treatment[temp$agency == "dod"] <-
    as.character(qualtrics[qualtrics$ID == IDs[i], "exp_1_gender_dod2"])
  temp$gender_treatment[temp$agency == "treasury"] <-
    as.character(qualtrics[qualtrics$ID == IDs[i], "exp_1_gender_treas2"])
  temp$gender_treatment[temp$agency == "hhs"] <-
    as.character(qualtrics[qualtrics$ID == IDs[i], "exp_1_gender_hhs2"])
  temp$gender_treatment[temp$agency == "educ"] <-
    as.character(qualtrics[qualtrics$ID == IDs[i], "exp_1_gender_educ2"])

  temp$dose[temp$agency == "dod"] <-
    as.character(qualtrics[qualtrics$ID == IDs[i], "pct_dose_dod"])
  temp$dose[temp$agency == "treasury"] <-
    as.character(qualtrics[qualtrics$ID == IDs[i], "pct_dose_treas"])
  temp$dose[temp$agency == "hhs"] <-
    as.character(qualtrics[qualtrics$ID == IDs[i], "pct_dose_hhs"])
  temp$dose[temp$agency == "educ"] <-
    as.character(qualtrics[qualtrics$ID == IDs[i], "pct_dose_educ"])

  # outcomes
  temp$ability[temp$agency == "dod"] <-
    qualtrics[qualtrics$ID == IDs[i], "dod_ability_dept2"]
  temp$ability[temp$agency == "treasury"] <-
    qualtrics[qualtrics$ID == IDs[i], "treasury_abil_dep2"]
  temp$ability[temp$agency == "hhs"] <-
    qualtrics[qualtrics$ID == IDs[i], "hhs_ability_dept2"]
  temp$ability[temp$agency == "educ"] <-
    qualtrics[qualtrics$ID == IDs[i], "educ_ability_dept2"]

  temp$rep_women[temp$agency == "dod"] <-
    qualtrics[qualtrics$ID == IDs[i], "dod_rep_women2"]
  temp$rep_women[temp$agency == "treasury"] <-
    qualtrics[qualtrics$ID == IDs[i], "treasury_rep_women2"]
  temp$rep_women[temp$agency == "hhs"] <-
    qualtrics[qualtrics$ID == IDs[i], "hhs_rep_women2"]
  temp$rep_women[temp$agency == "educ"] <-
    qualtrics[qualtrics$ID == IDs[i], "educ_rep_women2"]

  if (qualtrics[qualtrics$ID == IDs[i], "ag_order1"] == "dod") {
    temp$ag_order[temp$agency == "dod"] <- 1
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order1"] == "hhs") {
    temp$ag_order[temp$agency == "hhs"] <- 1
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order1"] == "educ") {
    temp$ag_order[temp$agency == "educ"] <- 1
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order1"] == "treasury") {
    temp$ag_order[temp$agency == "treasury"] <- 1
  }

  if (qualtrics[qualtrics$ID == IDs[i], "ag_order2"] == "dod") {
    temp$ag_order[temp$agency == "dod"] <- 2
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order2"] == "hhs") {
    temp$ag_order[temp$agency == "hhs"] <- 2
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order2"] == "educ") {
    temp$ag_order[temp$agency == "educ"] <- 2
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order2"] == "treasury") {
    temp$ag_order[temp$agency == "treasury"] <- 2
  }

  if (qualtrics[qualtrics$ID == IDs[i], "ag_order3"] == "dod") {
    temp$ag_order[temp$agency == "dod"] <- 3
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order3"] == "hhs") {
    temp$ag_order[temp$agency == "hhs"] <- 3
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order3"] == "educ") {
    temp$ag_order[temp$agency == "educ"] <- 3
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order3"] == "treasury") {
    temp$ag_order[temp$agency == "treasury"] <- 3
  }

  if (qualtrics[qualtrics$ID == IDs[i], "ag_order4"] == "dod") {
    temp$ag_order[temp$agency == "dod"] <- 4
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order4"] == "hhs") {
    temp$ag_order[temp$agency == "hhs"] <- 4
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order4"] == "educ") {
    temp$ag_order[temp$agency == "educ"] <- 4
  }
  if (qualtrics[qualtrics$ID == IDs[i], "ag_order4"] == "treasury") {
    temp$ag_order[temp$agency == "treasury"] <- 4
  }


  temp0 <- rbind.data.frame(temp0, temp)
}
head(temp0)

# remove empty rows
dim(temp0)
qualtrics1 <- temp0[temp0$ID != 0, ]
head(qualtrics1)
dim(qualtrics1)

# merge in additional respondent-level variables. Note "rank" treatment same for all four questions/independently randomized.
dim(qualtrics1)
qualtrics1 <- merge(qualtrics1, qualtrics[, c(
  "ID",
  "female",
  "party2",
  "rank",
  "nhwhite",
  "nhblack",
  "hispanic",
  "nhasian",
  "ideo2",
  "order_exp"
)], by = "ID", all.x = T)
dim(qualtrics1)
table(qualtrics1$rank)
head(qualtrics1)

class(qualtrics1$party2)
qualtrics1$party2 <- factor(qualtrics1$party2, levels = c("rep", "ind", "dem"))
table(qualtrics1$party2)

class(qualtrics1$dose)
qualtrics1$dose <-
  factor(qualtrics1$dose,
    levels = c("control", "20", "30", "40", "50", "60", "70", "80")
  )
table(qualtrics1$dose)

table(qualtrics1$gender_treatment, exclude = NULL)
qualtrics1$gender_treatment[qualtrics1$dose == "control"] <- "control"
qualtrics1$gender_treatment <-
  factor(qualtrics1$gender_treatment, levels = c("men", "control", "women"))
table(qualtrics1$gender_treatment)

head(qualtrics1)

# code dose treatments in terms of implied pairs (%women, %men)
qualtrics1$dose2 <- NA
qualtrics1$dose2[qualtrics1$dose == "20" & qualtrics1$gender_treatment == "women"] <- "20_80"
qualtrics1$dose2[qualtrics1$dose == "80" & qualtrics1$gender_treatment == "men"] <- "20_80"

qualtrics1$dose2[qualtrics1$dose == "30" & qualtrics1$gender_treatment == "women"] <- "30_70"
qualtrics1$dose2[qualtrics1$dose == "70" & qualtrics1$gender_treatment == "men"] <- "30_70"

qualtrics1$dose2[qualtrics1$dose == "40" & qualtrics1$gender_treatment == "women"] <- "40_60"
qualtrics1$dose2[qualtrics1$dose == "60" & qualtrics1$gender_treatment == "men"] <- "40_60"

qualtrics1$dose2[qualtrics1$dose == "50" & qualtrics1$gender_treatment == "women"] <- "50_50"
qualtrics1$dose2[qualtrics1$dose == "50" & qualtrics1$gender_treatment == "men"] <- "50_50"

qualtrics1$dose2[qualtrics1$dose == "60" & qualtrics1$gender_treatment == "women"] <- "60_40"
qualtrics1$dose2[qualtrics1$dose == "40" & qualtrics1$gender_treatment == "men"] <- "60_40"

qualtrics1$dose2[qualtrics1$dose == "70" & qualtrics1$gender_treatment == "women"] <- "70_30"
qualtrics1$dose2[qualtrics1$dose == "30" & qualtrics1$gender_treatment == "men"] <- "70_30"

qualtrics1$dose2[qualtrics1$dose == "80" & qualtrics1$gender_treatment == "women"] <- "80_20"
qualtrics1$dose2[qualtrics1$dose == "20" & qualtrics1$gender_treatment == "men"] <- "80_20"

qualtrics1$dose2[qualtrics1$dose == "control"] <- "control"
table(qualtrics1$dose2, exclude = NULL)

head(qualtrics1)
class(qualtrics1$ability)
class(qualtrics1$rep_women)



## spotcheck that long data are assembled correctly
check <- c(T)
for (i in 1:1000) {
  temp <- sample(IDs, size = 1)
  check <- c(check, qualtrics1$rep_women[qualtrics1$agency == "hhs" &
    qualtrics1$ID == IDs[i]] == qualtrics$hhs_rep_women2[qualtrics$ID ==
    IDs[i]])
  check <- c(check, qualtrics1$rep_women[qualtrics1$agency == "dod" &
    qualtrics1$ID == IDs[i]] == qualtrics$dod_rep_women2[qualtrics$ID ==
    IDs[i]])
  check <- c(check, qualtrics1$rep_women[qualtrics1$agency == "educ" &
    qualtrics1$ID == IDs[i]] == qualtrics$educ_rep_women2[qualtrics$ID ==
    IDs[i]])
  check <- c(check, qualtrics1$rep_women[qualtrics1$agency == "treasury" &
    qualtrics1$ID == IDs[i]] == qualtrics$treasury_rep_women2[qualtrics$ID ==
    IDs[i]])

  check <- c(check, qualtrics1$ability[qualtrics1$agency == "hhs" &
    qualtrics1$ID == IDs[i]] == qualtrics$hhs_ability2[qualtrics$ID == IDs[i]])
  check <- c(check, qualtrics1$ability[qualtrics1$agency == "dod" &
    qualtrics1$ID == IDs[i]] == qualtrics$dod_ability2[qualtrics$ID == IDs[i]])
  check <- c(check, qualtrics1$ability[qualtrics1$agency == "educ" &
    qualtrics1$ID == IDs[i]] == qualtrics$educ_ability2[qualtrics$ID == IDs[i]])
  check <- c(check, qualtrics1$ability[qualtrics1$agency == "treasury" &
    qualtrics1$ID == IDs[i]] == qualtrics$treasury_ability2[qualtrics$ID ==
    IDs[i]])
}
table(check)




# Descriptive statistics -------------------------


### Table B2 ----------------------------------------

sum.ability <- c(summary(qualtrics1$ability), sd(qualtrics1$ability))
sum.rep_women <- c(summary(qualtrics1$rep_women), sd(qualtrics1$rep_women))
names(sum.rep_women) <-
  c(names(sum.rep_women)[1:(length(names(sum.rep_women)) - 1)], "Std. Dev.")
vars <- names(sum.rep_women)
temp <- rbind.data.frame(sum.ability, sum.rep_women)
colnames(temp) <- vars
rownames(temp) <- c(
  "Agency Fulfill Mission",
  "Agency Represent Women"
)

temp.01 <- temp
temp.01 <- round(temp.01, 3)
temp.01

print(xtable(temp.01))

## put back on 1-4 scale for interpretation
temp[1, ] <- temp[1, ] * 4 - temp[1, ] + 1
temp[2, ] <- temp[2, ] * 4 - temp[2, ] + 1

temp <- round(temp, 3)
temp




# ANALYSIS ----------------------------------------------------



### Table 1 ------------------------------------------


class(qualtrics1$party2)
qualtrics1$party2 <- factor(qualtrics1$party2, levels = c("rep", "ind", "dem"))
levels(qualtrics1$party2)

## POOLED ANALYSIS FOR DOSE RESPONSE
# rep women pooled
m1 <- lm(rep_women ~ gender_treatment, data = qualtrics1)
temp1 <- coeftest(m1, vcov = vcovCL, cluster = ~ID)
temp1
m1$se <- temp1[, "Std. Error"]

# rep women with party interaction
m2 <- lm(rep_women ~ gender_treatment * party2, data = qualtrics1)
temp2 <- coeftest(m2, vcov = vcovCL, cluster = ~ID)
temp2
m2$se <- temp2[, "Std. Error"]

# rep women with gender interaction
m3 <- lm(rep_women ~ gender_treatment * female, data = qualtrics1)
temp3 <- coeftest(m3, vcov = vcovCL, cluster = ~ID)
temp3
m3$se <- temp3[, "Std. Error"]

## ability pooled
m4 <- lm(ability ~ gender_treatment, data = qualtrics1)
temp4 <- coeftest(m4, vcov = vcovCL, cluster = ~ID)
temp4
m4$se <- temp4[, "Std. Error"]
m4$coefficients
m4$se

## ability with party interaction
m5 <- lm(ability ~ gender_treatment * party2, data = qualtrics1)
temp5 <- coeftest(m5, vcov = vcovCL, cluster = ~ID)
temp5
m5$se <- temp5[, "Std. Error"]

## ability with gender interaction
m6 <- lm(ability ~ gender_treatment * female, data = qualtrics1)
temp6 <- coeftest(m6, vcov = vcovCL, cluster = ~ID)
temp6
m6$se <- temp6[, "Std. Error"]


stargazer(
  m1,
  m2,
  m3,
  m4,
  m5,
  m6,
  dep.var.labels = c("Represent Women", "Fulfill Mission"),
  title = "Effects of Gender Treatment in Dose Response Experiment. The table below shows
          average treatment effects of emphasizing women (relative to men) in the dose response experiment pooled
          across all doses. Models (a) and (d) show average effects in the entire sample.
          Models (b)-(c) and (e)-(f) condition on respondent party and gender, respectively.",
  covariate.labels = c(
    "condition: control",
    "condition: women",
    "Independent",
    "Democrat",
    "control * Indepdenent",
    "women * Independent",
    "control * Democrat",
    "women * Democrat",
    "female respondent",
    "control * female",
    "women * female",
    "(Intercept)"
  ),
  se = list(m1$se, m2$se, m3$se, m4$se, m5$se, m6$se),
  label = "table:dose_pooled",
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)




### Figure 3 ------------------------------------


# Mean response to Represent Women Question by Dose and Gender Treatment

# estimate means
res2 <- dose_means(
  data = qualtrics1,
  dv = "rep_women",
  color_women = "#D55E00",
  color_men = "#0072B2",
  shape_women = 19,
  shape_men = 17
)
head(res2)
control <- res2[res2$gender == "control", ] # extract control results
res2 <- res2[res2$gender != "control", ] # keep all other conditions
res2$order <- 1:nrow(res2)
res2 <- res2[rev(res2$order), ]
groups <- list(c(1, 2), c(3, 4), c(5, 6), c(7, 8), c(9, 10), c(11, 12), c(13, 14))
diffs <- NA
res2$weight <- NA

# test whether pooled regression returns same as weighted sum of treatment effects across bins
for (i in 1:length(groups)) {
  num1 <- groups[[i]][1]
  num2 <- groups[[i]][2]


  n1 <- nrow(qualtrics1[qualtrics1$gender_treatment == res2[num1, "gender"], ])
  n2 <- nrow(qualtrics1[qualtrics1$gender_treatment == res2[num2, "gender"], ])

  temp1 <- nrow(qualtrics1[qualtrics1$gender_treatment == res2[num1, "gender"] &
    qualtrics1$dose == res2[num1, "dose_test"] &
    qualtrics1$gender_treatment != "control", ])
  res2$weight[num1] <- temp1 / n1

  temp2 <- nrow(qualtrics1[qualtrics1$gender_treatment == res2[num2, "gender"] &
    qualtrics1$dose == res2[num2, "dose_test"] &
    qualtrics1$gender_treatment != "control", ])
  res2$weight[num2] <- temp2 / n2

  diffs[i] <-
    res2$coef[num1] * res2$weight[num1] - res2$coef[num2] * res2$weight[num2]
}

diffs
sum(diffs)
summary(lm(rep_women ~ gender_treatment, data = qualtrics1[qualtrics1$gender_treatment !=
  "control", ]))$coefficients["gender_treatmentwomen", "Estimate"]

## same result


# Figure 3
pdf(file = "Accepted Paper Replication Archive/figures/figure_3.pdf")
par(mar = c(6, 4, 4, 4))
plot(
  1:7,
  res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"],
  ylim = c(min(c(res2$lb, control$lb)) - .03, max(res2$ub)),
  axes = F,
  xlab = "",
  ylab = "Represent Women's Interests (0-1 Scale)",
  col = res2$color[res2$gender == "women"]
)
polygon(
  x = c(0, 8, 8, 0),
  y = c(control$lb, control$lb, control$ub, control$ub),
  col = t_col("orange", percent = 60),
  border = "white"
)
segments(
  x0 = 0,
  x1 = 8,
  y0 = control$coef,
  y1 = control$coef,
  col = "darkorange",
  lty = 2
)
text(
  x = 4,
  y = control$lb + .005,
  labels = "Control",
  col = "darkorange"
)
points(1:7, res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"], col = res2$color[res2$gender ==
    "women"]
)
segments(
  x0 = 1:7,
  x1 = 1:7,
  y0 = res2$lb[res2$gender == "women"],
  y1 = res2$ub[res2$gender == "women"],
  col = res2$color[res2$gender == "women"]
)
points(c(1:7) + .1,
  res2$coef[res2$gender == "men"],
  pch = res2$pch[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
segments(
  x0 = c(1:7) + .1,
  x1 = c(1:7) + .1,
  y0 = res2$lb[res2$gender == "men"],
  y1 = res2$ub[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
axis(
  1,
  at = c(1, 3, 5, 7),
  labels = rev(
    c(
      "20%Men\n80%Women",
      "40%Men\n60%Women",
      "60%Men\n40%Women",
      "80%Men\n20%Women"
    )
  ),
  cex.axis = .8,
  las = 1,
  tck = -.02
)
axis(1,
  at = c(2, 4, 6),
  tck = .02,
  labels = rep("", 3)
)

text(
  x = c(2, 4, 6),
  y = c(min(c(res2$lb, control$lb)) - .025),
  rev(
    c("30%Men\n70%Women", "50%Men\n50%Women", "70%Men\n30%Women")
  ),
  cex = .8
)

axis(2, at = seq(0, 1, .05), las = 2)

legend(
  "topleft",
  pch = c(19, 17),
  lty = 1,
  legend = c(
    "% women treatment",
    "% men treatment"
  ),
  col = c(
    res2$color[res2$gender == "women"][1],
    res2$color[res2$gender == "men"][1]
  ),
  cex = 1
)
dev.off()


### Figure 4 --------------------------------------------------------
# ABILIITY (run agency well)

res2 <- dose_means(
  data = qualtrics1,
  dv = "ability",
  color_women = "#D55E00",
  color_men = "#0072B2",
  shape_women = 19,
  shape_men = 17
)
res2
control <- res2[res2$gender == "control", ]
res2 <- res2[res2$gender != "control", ]
res2$order <- 1:nrow(res2)
res2 <- res2[rev(res2$order), ]


# Figure 4
pdf(file = "Accepted Paper Replication Archive/figures/figure_4.pdf")

par(mar = c(6, 4, 4, 4))
plot(
  1:7,
  res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"],
  ylim = c(min(c(res2$lb, control$lb)) - .03, max(res2$ub)),
  axes = F,
  xlab = "",
  ylab = "Agency Will Fulfill Mission (0-1 Scale)",
  col = res2$color[res2$gender == "women"]
)
polygon(
  x = c(0, 8, 8, 0),
  y = c(control$lb, control$lb, control$ub, control$ub),
  col = t_col("orange", percent = 60),
  border = "white"
)
segments(
  x0 = 0,
  x1 = 8,
  y0 = control$coef,
  y1 = control$coef,
  col = "darkorange",
  lty = 2
)
text(
  x = 4,
  y = control$lb + .005,
  labels = "Control",
  col = "darkorange"
)
points(1:7, res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"], col = res2$color[res2$gender ==
    "women"]
)
segments(
  x0 = 1:7,
  x1 = 1:7,
  y0 = res2$lb[res2$gender == "women"],
  y1 = res2$ub[res2$gender == "women"],
  col = res2$color[res2$gender == "women"]
)
points(c(1:7) + .1,
  res2$coef[res2$gender == "men"],
  pch = res2$pch[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
segments(
  x0 = c(1:7) + .1,
  x1 = c(1:7) + .1,
  y0 = res2$lb[res2$gender == "men"],
  y1 = res2$ub[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
axis(
  1,
  at = c(1, 3, 5, 7),
  labels = rev(
    c(
      "20%Men\n80%Women",
      "40%Men\n60%Women",
      "60%Men\n40%Women",
      "80%Men\n20%Women"
    )
  ),
  cex.axis = .8,
  las = 1,
  tck = -.02
)
axis(1,
  at = c(2, 4, 6),
  tck = .02,
  labels = rep("", 3)
)

text(
  x = c(2, 4, 6),
  y = c(min(c(res2$lb, control$lb)) - .025),
  rev(
    c("30%Men\n70%Women", "50%Men\n50%Women", "70%Men\n30%Women")
  ),
  cex = .8
)

axis(2, at = seq(0, 1, .05), las = 2)

legend(
  "topleft",
  pch = c(19, 17),
  lty = 1,
  legend = c(
    "% women treatment",
    "% men treatment"
  ),
  col = c(
    res2$color[res2$gender == "women"][1],
    res2$color[res2$gender == "men"][1]
  ),
  cex = 1
)
dev.off()






# Robustness Checks  -----------------------------------------

### Figures B10-B11  -----------------------------------------



## Represent Women
dose_test <-
  c("80_20", "70_30", "60_40", "50_50", "40_60", "30_70", "20_80")

res2 <- as.data.frame(matrix(nrow = 7, ncol = 3))
colnames(res2) <- c("dose", "coef", "se")
res2$dose <- dose_test
res2

for (i in 1:length(dose_test)) {
  m <-
    lm(rep_women ~ gender_treatment, data = qualtrics1[qualtrics1$dose2 == c(dose_test[i]), ])
  temp <- coeftest(m, vcov = vcovCL, cluster = ~ID)
  res2$coef[res2$dose == dose_test[i]] <- temp[2, "Estimate"]
  res2$se[res2$dose == dose_test[i]] <- temp[2, "Std. Error"]
}
res2$lb <- res2$coef - 1.96 * res2$se
res2$ub <- res2$coef + 1.96 * res2$se
res2
# tau_rep_women_qualtrics.pdf
pdf(file = "Accepted Paper Replication Archive/figures/figure_b10.pdf")
plot(
  1:7,
  rev(res2$coef),
  ylim = c(-.06, .1),
  pch = 19,
  col = "black",
  ylab = "Difference in Effect Between Women/Men Condition",
  axes = F,
  xlab = ""
)
segments(1:7, rev(res2$lb), 1:7, rev(res2$ub), lty = 1)
abline(h = 0, col = "red", lty = 2)
axis(
  1,
  at = c(1, 3, 5, 7),
  labels = rev(
    c(
      "20%Men\n80%Women",
      "40%Men\n60%Women",
      "60%Men\n40%Women",
      "80%Men\n20%Women"
    )
  ),
  cex.axis = .8,
  las = 1,
  tck = -.02
)
axis(1,
  at = c(2, 4, 6),
  tck = .02,
  labels = rep("", 3)
)

text(
  x = c(2, 4, 6),
  y = c(min(c(res2$lb)) - .045),
  rev(
    c("30%Men\n70%Women", "50%Men\n50%Women", "70%Men\n30%Women")
  ),
  cex = .8
)
axis(2, at = seq(-1, 1, by = .02), las = 2)
dev.off()

### run agency well
dose_test <-
  c("80_20", "70_30", "60_40", "50_50", "40_60", "30_70", "20_80")

res2 <- as.data.frame(matrix(nrow = 7, ncol = 3))
colnames(res2) <- c("dose", "coef", "se")
res2$dose <- dose_test
res2

for (i in 1:length(dose_test)) {
  m <-
    lm(ability ~ gender_treatment, data = qualtrics1[qualtrics1$dose2 == c(dose_test[i]), ])
  temp <- coeftest(m, vcov = vcovCL, cluster = ~ID)
  res2$coef[res2$dose == dose_test[i]] <- temp[2, "Estimate"]
  res2$se[res2$dose == dose_test[i]] <- temp[2, "Std. Error"]
}
res2$lb <- res2$coef - 1.96 * res2$se
res2$ub <- res2$coef + 1.96 * res2$se
res2

# Figure B11
pdf(file = "Accepted Paper Replication Archive/figures/figure_b11.pdf")

plot(
  1:7,
  rev(res2$coef),
  ylim = c(-.06, .1),
  pch = 19,
  col = "black",
  ylab = "Difference in Effect Between Women/Men Condition",
  axes = F,
  xlab = ""
)
segments(1:7, rev(res2$lb), 1:7, rev(res2$ub), lty = 1)
abline(h = 0, col = "red", lty = 2)
axis(
  1,
  at = c(1, 3, 5, 7),
  labels = rev(
    c(
      "20%Men\n80%Women",
      "40%Men\n60%Women",
      "60%Men\n40%Women",
      "80%Men\n20%Women"
    )
  ),
  cex.axis = .8,
  las = 1,
  tck = -.02
)
axis(1,
  at = c(2, 4, 6),
  tck = .02,
  labels = rep("", 3)
)

text(
  x = c(2, 4, 6),
  y = c(min(c(res2$lb)) - .035),
  rev(
    c("30%Men\n70%Women", "50%Men\n50%Women", "70%Men\n30%Women")
  ),
  cex = .8
)
axis(2, at = seq(-1, 1, by = .02), las = 2)

dev.off()



### Table B3  -----------------------------------------


### condition on rank of job described and agency
# represent women

m1 <- lm(rep_women ~ gender_treatment + rank + factor(agency), data = qualtrics1)
temp1 <- coeftest(m1, vcov = vcovCL, cluster = ~ID)
m1$se <- temp1[, "Std. Error"]

m2 <- lm(rep_women ~ gender_treatment + rank + gender_treatment * factor(agency), data = qualtrics1)
temp2 <- coeftest(m2, vcov = vcovCL, cluster = ~ID)
m2$se <- temp2[, "Std. Error"]

m3 <- lm(rep_women ~ gender_treatment + rank + gender_treatment * rank + factor(agency), data = qualtrics1)
temp3 <- coeftest(m3, vcov = vcovCL, cluster = ~ID)
m3$se <- temp3[, "Std. Error"]

# run agency well

m4 <- lm(ability ~ gender_treatment + rank + factor(agency), data = qualtrics1)
temp4 <- coeftest(m4, vcov = vcovCL, cluster = ~ID)
m4$se <- temp4[, "Std. Error"]

m5 <- lm(ability ~ gender_treatment + rank + gender_treatment * factor(agency), data = qualtrics1)
temp5 <- coeftest(m5, vcov = vcovCL, cluster = ~ID)
m5$se <- temp5[, "Std. Error"]

m6 <- lm(ability ~ gender_treatment + rank + gender_treatment * rank + factor(agency), data = qualtrics1)
temp6 <- coeftest(m6, vcov = vcovCL, cluster = ~ID)
m6$se <- temp6[, "Std. Error"]


stargazer(
  m1,
  m2,
  m3,
  m4,
  m5,
  m6,
  dep.var.labels = c("Represent Women", "Fulfill Mission"),
  title = "\\textbf{Condition on Agency and Rank.} Effects of Gender Treatment in Dose Response Experiment.",
  se = list(m1$se, m2$se, m3$se, m4$se, m5$se, m6$se),
  label = "table:dose_pooled_agency_rank",
  omit.stat = c("f", "rsq", "adj.rsq", "ser"),
  covariate.labels = c(
    "condition: control",
    "condition: women",
    "Rank: `top jobs'",
    "Agency: Education",
    "Agency: HHS",
    "Agency: Treasury",
    "control * Education",
    "women * Education",
    "control * HHS",
    "women * HHS",
    "control * Treasury",
    "women * Treasury",
    "control * Rank: `top jobs'",
    "women * Rank: `top jobs'"
  )
)



### Table B4  -----------------------------------------

# represent women, respondent gender
models <- list(NA)
ses <- list(NA)
dose_test
for (i in 1:length(dose_test)) {
  temp0 <- lm(rep_women ~ gender_treatment * female, data = qualtrics1[qualtrics1$dose2 == dose_test[i], ])
  temp1 <- coeftest(temp0, vcov = vcovCL, cluster = ~ID)
  temp0$se <- temp1[, "Std. Error"]

  ses[[i]] <- temp0$se
  models[[i]] <- temp0
}
names(ses) <- dose_test
names(models) <- dose_test



stargazer(models,
  dep.var.labels = "Represent Women",
  column.labels = c("80/20", "70/30", "60/40", "50/50", "40/60", "30/70", "20/80"),
  title = "\\textbf{Treatment Effects by Dose, Interacted with Respondent Gender: Agency will Represent Women's Interests}.
  Column labels indicate percent men / percent women among respondents in each model. Coefficient on `women' condition indicates effect of providing information on percent women relative to a logically equivalent message about percent men.",
  se = ses,
  label = "table:dose_female_separate_rep_women",
  covariate.labels = c(
    "condition: women",
    "female respondent",
    "women * female respondent",
    "(Intercept)"
  ),
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)


### Table B5  -----------------------------------------


## run agency well, respondent gender
models <- list(NA)
ses <- list(NA)
dose_test
for (i in 1:length(dose_test)) {
  temp0 <- lm(ability ~ gender_treatment * female, data = qualtrics1[qualtrics1$dose2 == dose_test[i], ])
  temp1 <- coeftest(temp0, vcov = vcovCL, cluster = ~ID)
  temp0$se <- temp1[, "Std. Error"]

  ses[[i]] <- temp0$se
  models[[i]] <- temp0
}
names(ses) <- dose_test
names(models) <- dose_test



stargazer(models,
  dep.var.labels = "Agency Fulfill Mission",
  column.labels = c("80/20", "70/30", "60/40", "50/50", "40/60", "30/70", "20/80"),
  title = "\\textbf{Treatment Effects by Dose, Interacted with Respondent Gender: Agency will Fulfill Mission}.
  Column labels indicate percent men / percent women among respondents in each model. Coefficient on `women' condition indicates effect of providing information on percent women relative to a logically equivalent message about percent men.",
  se = ses,
  label = "table:dose_female_separate_ability",
  covariate.labels = c(
    "condition: women",
    "female respondent",
    "women * female respondent",
    "(Intercept)"
  ),
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)


### Table B6  -----------------------------------------


## represent women, party
models <- list(NA)
ses <- list(NA)
dose_test
for (i in 1:length(dose_test)) {
  temp0 <- lm(rep_women ~ gender_treatment * party2, data = qualtrics1[qualtrics1$dose2 == dose_test[i], ])
  temp1 <- coeftest(temp0, vcov = vcovCL, cluster = ~ID)
  temp0$se <- temp1[, "Std. Error"]

  ses[[i]] <- temp0$se
  models[[i]] <- temp0
}
names(ses) <- dose_test
names(models) <- dose_test



stargazer(models,
  dep.var.labels = "Represent Women",
  column.labels = c("80/20", "70/30", "60/40", "50/50", "40/60", "30/70", "20/80"),
  title = "\\textbf{Treatment Effects by Dose, Interacted with Respondent Party: Agency will Represent Women's Interests}.
  Column labels indicate percent men / percent women among respondents in each model. Coefficient on `women' condition indicates effect of providing information on percent women relative to a logically equivalent message about percent men.",
  se = ses,
  label = "table:dose_party_separate_rep_women",
  covariate.labels = c(
    "condition: women",
    "Independent",
    "Democrat",
    "women * Independent",
    "women * Democrat",
    "(Intercept)"
  ),
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)


### Table B7  -----------------------------------------


## run agency well, party
models <- list(NA)
ses <- list(NA)
dose_test
for (i in 1:length(dose_test)) {
  temp0 <- lm(ability ~ gender_treatment * party2, data = qualtrics1[qualtrics1$dose2 == dose_test[i], ])
  temp1 <- coeftest(temp0, vcov = vcovCL, cluster = ~ID)
  temp0$se <- temp1[, "Std. Error"]

  ses[[i]] <- temp0$se
  models[[i]] <- temp0
}
names(ses) <- dose_test
names(models) <- dose_test



stargazer(models,
  dep.var.labels = "Agency Fulfill Mission",
  column.labels = c("80/20", "70/30", "60/40", "50/50", "40/60", "30/70", "20/80"),
  title = "\\textbf{Treatment Effects by Dose, Interacted with Respondent Party: Agency will Fulfill Mission}.
  Column labels indicate percent men / percent women among respondents in each model. Coefficient on `women' condition indicates effect of providing information on percent women relative to a logically equivalent message about percent men.",
  se = ses,
  label = "table:dose_party_separate_ability",
  covariate.labels = c(
    "condition: women",
    "Independent",
    "Democrat",
    "women * Independent",
    "women * Democrat",
    "(Intercept)"
  ),
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)




### Order Effects



### Table B8  -----------------------------------------
## Condition on order in which experiments were presented

table(qualtrics1$order_exp)
qualtrics1$order_exp2 <- factor(qualtrics1$order_exp, levels = c("job", "dose", "charts"))
qualtrics1$order_exp3 <- factor(qualtrics1$order_exp, levels = c("charts", "dose", "job"))

m1 <- lm(rep_women ~ gender_treatment * order_exp, data = qualtrics1)
temp1 <- coeftest(m1, vcov = vcovCL, cluster = ~ID)
m1$se <- temp1[, "Std. Error"]

m2 <- lm(ability ~ gender_treatment * order_exp, data = qualtrics1)
temp2 <- coeftest(m2, vcov = vcovCL, cluster = ~ID)
m2$se <- temp2[, "Std. Error"]



stargazer(m1, m2,
  column.labels = c("Represent Women", "Agency Fulfill Mission"),
  dep.var.labels = c("", ""),
  covariate.labels = c("condition: control", "condition: women", "nominee experiment first", "visualization experiment first", "control*nominee first", "women*nominee first", "control*visualization first", "women*visualization first"),
  title = "\\textbf{Effect of Experiment Order on Treatment Effects.} Omitted category is condition presenting framing experiment first in the survey instrument.",
  omit.stat = c("f", "rsq", "adj.rsq", "ser"),
  se = list(m1$se, m2$se),
  label = "app:table_order_effects"
)





### Figures B12-B14  -----------------------------------------
### treatment effects conditional on which experiment presented first



## Figure B12
## rep women dose first
res2 <- dose_means(
  data = qualtrics1[qualtrics1$order_exp == "dose", ], dv = "rep_women", color_women = "#D55E00",
  color_men = "#0072B2", shape_women = 19, shape_men = 17
)
head(res2)
control <- res2[res2$gender == "control", ]
res2 <- res2[res2$gender != "control", ]
res2$order <- 1:nrow(res2)
res2 <- res2[rev(res2$order), ]

# a_dose_rep_women_qualtrics_2_dosefirst.pdf
pdf(file = "Accepted Paper Replication Archive/figures/figure_b12.pdf")

par(mar = c(6, 4, 4, 4))
plot(1:7, res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"], ylim = c(min(c(res2$lb, control$lb)) - .03, max(res2$ub)), axes = F,
  xlab = "", ylab = "Represent Women's Interests (0-1 Scale)",
  col = res2$color[res2$gender == "women"], main = "Framing Experiment Presented First"
)
polygon(
  x = c(0, 8, 8, 0), y = c(control$lb, control$lb, control$ub, control$ub),
  col = t_col("orange", percent = 60), border = "white"
)
segments(x0 = 0, x1 = 8, y0 = control$coef, y1 = control$coef, col = "darkorange", lty = 2)
text(x = 4, y = control$lb + .005, labels = "Control", col = "darkorange")
points(1:7, res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"], col = res2$color[res2$gender == "women"]
)
segments(
  x0 = 1:7, x1 = 1:7, y0 = res2$lb[res2$gender == "women"],
  y1 = res2$ub[res2$gender == "women"],
  col = res2$color[res2$gender == "women"]
)
points(c(1:7) + .1, res2$coef[res2$gender == "men"],
  pch = res2$pch[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
segments(
  x0 = c(1:7) + .1, x1 = c(1:7) + .1, y0 = res2$lb[res2$gender == "men"],
  y1 = res2$ub[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
axis(1, at = c(1, 3, 5, 7), labels = rev(c(
  "20%Men\n80%Women", "40%Men\n60%Women",
  "60%Men\n40%Women",
  "80%Men\n20%Women"
)), cex.axis = .8, las = 1, tck = -.02)
axis(1, at = c(2, 4, 6), tck = .02, labels = rep("", 3))

text(
  x = c(2, 4, 6), y = c(min(c(res2$lb, control$lb)) - .025),
  rev(c("30%Men\n70%Women", "50%Men\n50%Women", "70%Men\n30%Women")),
  cex = .8
)

axis(2, at = seq(0, 1, .05), las = 2)


legend("topleft",
  pch = c(19, 17), lty = 1, legend = c(
    "% women treatment",
    "% men treatment"
  ),
  col = c(
    res2$color[res2$gender == "women"][1],
    res2$color[res2$gender == "men"][1]
  ), cex = 1
)
dev.off()


## Figure B13
## rep women, nominee experiment presented  first

res2 <- dose_means(
  data = qualtrics1[qualtrics1$order_exp == "job", ], dv = "rep_women", color_women = "#D55E00",
  color_men = "#0072B2", shape_women = 19, shape_men = 17
)
head(res2)
control <- res2[res2$gender == "control", ]
res2 <- res2[res2$gender != "control", ]
res2$order <- 1:nrow(res2)
res2 <- res2[rev(res2$order), ]

# a_dose_rep_women_qualtrics_2_jobfirst.pdf
pdf(file = "Accepted Paper Replication Archive/figures/figure_b13.pdf")

par(mar = c(6, 4, 4, 4))
plot(1:7, res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"], ylim = c(min(c(res2$lb, control$lb)) - .03, max(res2$ub)), axes = F,
  xlab = "", ylab = "Represent Women's Interests (0-1 Scale)",
  col = res2$color[res2$gender == "women"], main = "Nominee Experiment Presented First"
)
polygon(
  x = c(0, 8, 8, 0), y = c(control$lb, control$lb, control$ub, control$ub),
  col = t_col("orange", percent = 60), border = "white"
)
segments(x0 = 0, x1 = 8, y0 = control$coef, y1 = control$coef, col = "darkorange", lty = 2)
text(x = 4, y = control$lb + .005, labels = "Control", col = "darkorange")
points(1:7, res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"], col = res2$color[res2$gender == "women"]
)
segments(
  x0 = 1:7, x1 = 1:7, y0 = res2$lb[res2$gender == "women"],
  y1 = res2$ub[res2$gender == "women"],
  col = res2$color[res2$gender == "women"]
)
points(c(1:7) + .1, res2$coef[res2$gender == "men"],
  pch = res2$pch[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
segments(
  x0 = c(1:7) + .1, x1 = c(1:7) + .1, y0 = res2$lb[res2$gender == "men"],
  y1 = res2$ub[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
axis(1, at = c(1, 3, 5, 7), labels = rev(c(
  "20%Men\n80%Women", "40%Men\n60%Women",
  "60%Men\n40%Women",
  "80%Men\n20%Women"
)), cex.axis = .8, las = 1, tck = -.02)
axis(1, at = c(2, 4, 6), tck = .02, labels = rep("", 3))

text(
  x = c(2, 4, 6), y = c(min(c(res2$lb, control$lb)) - .025),
  rev(c("30%Men\n70%Women", "50%Men\n50%Women", "70%Men\n30%Women")),
  cex = .8
)

axis(2, at = seq(0, 1, .05), las = 2)


legend("topleft",
  pch = c(19, 17), lty = 1, legend = c(
    "% women treatment",
    "% men treatment"
  ),
  col = c(
    res2$color[res2$gender == "women"][1],
    res2$color[res2$gender == "men"][1]
  ), cex = 1
)
dev.off()




### Figure B14  -----------------------------------------
##### rep women charts first

res2 <- dose_means(
  data = qualtrics1[qualtrics1$order_exp == "charts", ], dv = "rep_women", color_women = "#D55E00",
  color_men = "#0072B2", shape_women = 19, shape_men = 17
)
head(res2)
control <- res2[res2$gender == "control", ]
res2 <- res2[res2$gender != "control", ]
res2$order <- 1:nrow(res2)
res2 <- res2[rev(res2$order), ]

# a_dose_rep_women_qualtrics_2_chartsfirst.pdf
pdf(file = "Accepted Paper Replication Archive/figures/figure_b14.pdf")

par(mar = c(6, 4, 4, 4))
plot(1:7, res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"], ylim = c(min(c(res2$lb, control$lb)) - .03, max(res2$ub)), axes = F,
  xlab = "", ylab = "Represent Women's Interests (0-1 Scale)",
  col = res2$color[res2$gender == "women"], main = "Visualization Experiment Presented First"
)
polygon(
  x = c(0, 8, 8, 0), y = c(control$lb, control$lb, control$ub, control$ub),
  col = t_col("orange", percent = 60), border = "white"
)
segments(x0 = 0, x1 = 8, y0 = control$coef, y1 = control$coef, col = "darkorange", lty = 2)
text(x = 4, y = control$lb + .005, labels = "Control", col = "darkorange")
points(1:7, res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"], col = res2$color[res2$gender == "women"]
)
segments(
  x0 = 1:7, x1 = 1:7, y0 = res2$lb[res2$gender == "women"],
  y1 = res2$ub[res2$gender == "women"],
  col = res2$color[res2$gender == "women"]
)
points(c(1:7) + .1, res2$coef[res2$gender == "men"],
  pch = res2$pch[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
segments(
  x0 = c(1:7) + .1, x1 = c(1:7) + .1, y0 = res2$lb[res2$gender == "men"],
  y1 = res2$ub[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
axis(1, at = c(1, 3, 5, 7), labels = rev(c(
  "20%Men\n80%Women", "40%Men\n60%Women",
  "60%Men\n40%Women",
  "80%Men\n20%Women"
)), cex.axis = .8, las = 1, tck = -.02)
axis(1, at = c(2, 4, 6), tck = .02, labels = rep("", 3))

text(
  x = c(2, 4, 6), y = c(min(c(res2$lb, control$lb)) - .025),
  rev(c("30%Men\n70%Women", "50%Men\n50%Women", "70%Men\n30%Women")),
  cex = .8
)

axis(2, at = seq(0, 1, .05), las = 2)


legend("topleft",
  pch = c(19, 17), lty = 1, legend = c(
    "% women treatment",
    "% men treatment"
  ),
  col = c(
    res2$color[res2$gender == "women"][1],
    res2$color[res2$gender == "men"][1]
  ), cex = 1
)
dev.off()


### Table B10  -----------------------------------------


## condition on order in which agency presented
# represent women

m1 <- lm(rep_women ~ gender_treatment * factor(ag_order), data = qualtrics1)
temp1 <- coeftest(m1, vcov = vcovCL, cluster = ~ID)
m1$se <- temp1[, "Std. Error"]

# run agency well
m2 <- lm(ability ~ gender_treatment * factor(ag_order), data = qualtrics1)
temp2 <- coeftest(m2, vcov = vcovCL, cluster = ~ID)
m2$se <- temp2[, "Std. Error"]

stargazer(
  m1,
  m2,
  dep.var.labels = c("Represent Women", "Fulfill Mission"),
  title = "\\textbf{Condition on order in which agencies were presented.} Effects of gender treatment in dose response experiment.",
  se = list(m1$se, m2$se),
  label = "table:dose_pooled_order_analysis",
  covariate.labels = c(
    "condition: control",
    "condition: women",
    "Order: 2",
    "Order: 3",
    "Order: 4",
    "control * Order: 2",
    "women * Order: 2",
    "control * Order: 3",
    "women * Order: 3",
    "control * Order: 4",
    "women * Order: 4",
    "(Intercept)"
  ),
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)


###
### Table B11  -----------------------------------------
###

## Limit to cases where dose response experiment presented first, and respondents are answering first question of four agencies presented.

## rep_women first response only and experiment comes first
head(qualtrics1)
m1_first_exp <-
  lm(rep_women ~ gender_treatment, data = qualtrics1[qualtrics1$ag_order == 1 &
    qualtrics1$order_exp == "dose", ])
temp1_first_exp <-
  coeftest(m1_first_exp, vcov = vcovCL, cluster = ~ID)
temp1_first_exp
m1_first_exp$se <- temp1_first_exp[, "Std. Error"]
m1_first_exp$coefficients
m1_first_exp$se

## ability first response only and experiment comes first
head(qualtrics1)
m2_first_exp <-
  lm(ability ~ gender_treatment, data = qualtrics1[qualtrics1$ag_order == 1 &
    qualtrics1$order_exp == "dose", ])
temp2_first_exp <-
  coeftest(m2_first_exp, vcov = vcovCL, cluster = ~ID)
temp2_first_exp
m2_first_exp$se <- temp2_first_exp[, "Std. Error"]
m2_first_exp$coefficients
m2_first_exp$se

### first response robustness table
stargazer(
  m1_first_exp,
  m2_first_exp,
  dep.var.labels = c("Represent Women", "Fulfill Mission"),
  title = "\\textbf{Robustness Check for Order Effects.} Effects of Gender Treatment in Dose Response Experiment. Estimates limited to the first of four items viewed by respondents, among respondents randomly assigned to participate in dose response experiment first in the survey instrument.",
  covariate.labels = c("condition: control", "condition: women", "(Intercept)"),
  se = list(m1_first_exp$se, m2_first_exp$se),
  label = "table:app_first_response_dose",
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)



##### Analysis of additional experiments in Qualtrics sample

# Nominee experiment  -----------------------------------------

## Exp 2
table(qualtrics$exp2_agency) ## agency treatment
qualtrics$e2_agency <- NA
qualtrics$e2_agency[qualtrics$exp2_agency == "Department of Health and Human Services"] <- "hhs"
qualtrics$e2_agency[qualtrics$exp2_agency == "Department of Defense"] <- "defense"
qualtrics$e2_agency[qualtrics$exp2_agency == "Department of Education"] <- "education"
qualtrics$e2_agency[qualtrics$exp2_agency == "Department of Treasury"] <- "treasury"
table(qualtrics$e2_agency, exclude = NULL)

table(qualtrics$exp2_gender) ## gender treatment
qualtrics$e2_gender <- NA
qualtrics$e2_gender[qualtrics$exp2_gender == "man"] <- "male"
qualtrics$e2_gender[qualtrics$exp2_gender == "woman"] <- "female"
qualtrics$e2_gender <- factor(qualtrics$e2_gender, levels = c("male", "female"))
table(qualtrics$e2_gender, exclude = NULL)



# Outcomes nominee experiment
## ability of nominee to effectively lead agency
table(qualtrics$exp2_dv1, exclude = NULL)
qualtrics$e2_nominee_confidence <- NA
qualtrics$e2_nominee_confidence[qualtrics$exp2_dv1 == "No confidence at all"] <- 1
qualtrics$e2_nominee_confidence[qualtrics$exp2_dv1 == "Not too much confidence"] <- 2
qualtrics$e2_nominee_confidence[qualtrics$exp2_dv1 == "Some confidence"] <- 3
qualtrics$e2_nominee_confidence[qualtrics$exp2_dv1 == "A lot of confidence"] <- 4
table(qualtrics$e2_nominee_confidence, exclude = NULL)

## president's ability to effectively staff government
table(qualtrics$exp2_dv2)
qualtrics$e2_pres_confidence <- NA
qualtrics$e2_pres_confidence[qualtrics$exp2_dv2 == "No confidence at all"] <- 1
qualtrics$e2_pres_confidence[qualtrics$exp2_dv2 == "Not too much confidence"] <- 2
qualtrics$e2_pres_confidence[qualtrics$exp2_dv2 == "Some confidence"] <- 3
qualtrics$e2_pres_confidence[qualtrics$exp2_dv2 == "A lot of confidence"] <- 4
table(qualtrics$e2_pres_confidence, exclude = NULL)

qualtrics$e2_nominee_confidence <- recode_0_1(qualtrics$e2_nominee_confidence)
qualtrics$e2_pres_confidence <- recode_0_1(qualtrics$e2_pres_confidence)



#  Chart Experiment  -----------------------------------------
table(qualtrics$FL_174_DO) ### chart displayed
qualtrics$e3_chart <- NA
qualtrics$e3_chart[qualtrics$FL_174_DO == "chart1"] <- "top_jobs_levels"
qualtrics$e3_chart[qualtrics$FL_174_DO == "chart2b"] <- "top_jobs_change"
qualtrics$e3_chart[qualtrics$FL_174_DO == "chart3"] <- "top_rf_levels"
qualtrics$e3_chart[qualtrics$FL_174_DO == "chart4"] <- "top_rf_US_levels"
qualtrics$e3_chart <- factor(qualtrics$e3_chart, levels = c(
  "top_jobs_change", "top_jobs_levels",
  "top_rf_levels",
  "top_rf_US_levels"
))
table(qualtrics$e3_chart, exclude = NULL)

## outcomes
## women have adequate rep in government
table(qualtrics$exp3_dv1)
table(qualtrics$exp3_dv1b)
table(qualtrics$exp3_dv1c)
table(qualtrics$exp3_dv1d)

table(qualtrics$exp3_dv3)
table(qualtrics$exp3_dv3b)
table(qualtrics$exp3_dv3c)
table(qualtrics$exp3_dv3d)

qualtrics$e3_women_rep <- NA
qualtrics$e3_women_rep[qualtrics$exp3_dv1 == "Much too little representation"] <- 1
qualtrics$e3_women_rep[qualtrics$exp3_dv1b == "Much too little representation"] <- 1
qualtrics$e3_women_rep[qualtrics$exp3_dv1c == "Much too little representation"] <- 1
qualtrics$e3_women_rep[qualtrics$exp3_dv1d == "Much too little representation"] <- 1

qualtrics$e3_women_rep[qualtrics$exp3_dv1 == "Too little representation"] <- 2
qualtrics$e3_women_rep[qualtrics$exp3_dv1b == "Too little representation"] <- 2
qualtrics$e3_women_rep[qualtrics$exp3_dv1c == "Too little representation"] <- 2
qualtrics$e3_women_rep[qualtrics$exp3_dv1d == "Too little representation"] <- 2

qualtrics$e3_women_rep[qualtrics$exp3_dv1 == "About the right amount of representation"] <- 3
qualtrics$e3_women_rep[qualtrics$exp3_dv1b == "About the right amount of representation"] <- 3
qualtrics$e3_women_rep[qualtrics$exp3_dv1c == "About the right amount of representation"] <- 3
qualtrics$e3_women_rep[qualtrics$exp3_dv1d == "About the right amount of representation"] <- 3

qualtrics$e3_women_rep[qualtrics$exp3_dv1 == "Too much representation"] <- 4
qualtrics$e3_women_rep[qualtrics$exp3_dv1b == "Too much representation"] <- 4
qualtrics$e3_women_rep[qualtrics$exp3_dv1c == "Too much representation"] <- 4
qualtrics$e3_women_rep[qualtrics$exp3_dv1d == "Too much representation"] <- 4

qualtrics$e3_women_rep[qualtrics$exp3_dv1 == "Much too much representation"] <- 5
qualtrics$e3_women_rep[qualtrics$exp3_dv1b == "Much too much representation"] <- 5
qualtrics$e3_women_rep[qualtrics$exp3_dv1c == "Much too much representation"] <- 5
qualtrics$e3_women_rep[qualtrics$exp3_dv1d == "Much too much representation"] <- 5
table(qualtrics$e3_women_rep, exclude = NULL)
qualtrics$e3_women_rep <- recode_0_1(qualtrics$e3_women_rep)
table(qualtrics$e3_women_rep, exclude = NULL)


qualtrics$e3_serve_US <- NA
qualtrics$e3_serve_US[qualtrics$exp3_dv3 == "No confidence at all"] <- 1
qualtrics$e3_serve_US[qualtrics$exp3_dv3b == "No confidence at all"] <- 1
qualtrics$e3_serve_US[qualtrics$exp3_dv3c == "No confidence at all"] <- 1
qualtrics$e3_serve_US[qualtrics$exp3_dv3d == "No confidence at all"] <- 1
table(qualtrics$e3_serve_US, exclude = NULL)


qualtrics$e3_serve_US[qualtrics$exp3_dv3 == "Not too much confidence"] <- 2
qualtrics$e3_serve_US[qualtrics$exp3_dv3b == "Not too much confidence"] <- 2
qualtrics$e3_serve_US[qualtrics$exp3_dv3c == "Not too much confidence"] <- 2
qualtrics$e3_serve_US[qualtrics$exp3_dv3d == "Not too much confidence"] <- 2
table(qualtrics$e3_serve_US, exclude = NULL)


qualtrics$e3_serve_US[qualtrics$exp3_dv3 == "Some confidence"] <- 3
qualtrics$e3_serve_US[qualtrics$exp3_dv3b == "Some confidence"] <- 3
qualtrics$e3_serve_US[qualtrics$exp3_dv3c == "Some confidence"] <- 3
qualtrics$e3_serve_US[qualtrics$exp3_dv3d == "Some confidence"] <- 3
table(qualtrics$e3_serve_US, exclude = NULL)


qualtrics$e3_serve_US[qualtrics$exp3_dv3 == "A lot of confidence"] <- 4
qualtrics$e3_serve_US[qualtrics$exp3_dv3b == "A lot of confidence"] <- 4
qualtrics$e3_serve_US[qualtrics$exp3_dv3c == "A lot of confidence"] <- 4
qualtrics$e3_serve_US[qualtrics$exp3_dv3d == "A lot of confidence"] <- 4
table(qualtrics$e3_serve_US, exclude = NULL)


table(qualtrics$e3_serve_US, exclude = NULL)
qualtrics$e3_serve_US <- recode_0_1(qualtrics$e3_serve_US)
table(qualtrics$e3_serve_US, exclude = NULL)


### Table B13  -----------------------------------------

class(qualtrics$party2)
qualtrics$party2 <- factor(qualtrics$party2, levels = c("rep", "ind", "dem"))
levels(qualtrics$party2)

m1 <- lm(e2_nominee_confidence ~ e2_gender, data = qualtrics)
temp <- coeftest(m1, vcov = vcovHC(m1, type = "HC1"))
temp
m1$se <- temp[, "Std. Error"]

m2 <- lm(e2_nominee_confidence ~ e2_gender * party2, data = qualtrics)
temp <- coeftest(m2, vcov = vcovHC(m2, type = "HC1"))
temp
m2$se <- temp[, "Std. Error"]

m3 <- lm(e2_nominee_confidence ~ e2_gender * female, data = qualtrics)
temp <- coeftest(m3, vcov = vcovHC(m3, type = "HC1"))
temp
m3$se <- temp[, "Std. Error"]

m4 <- lm(e2_pres_confidence ~ e2_gender, data = qualtrics)
temp <- coeftest(m4, vcov = vcovHC(m4, type = "HC1"))
temp
m4$se <- temp[, "Std. Error"]

m5 <- lm(e2_pres_confidence ~ e2_gender * party2, data = qualtrics)
temp <- coeftest(m5, vcov = vcovHC(m5, type = "HC1"))
temp
m5$se <- temp[, "Std. Error"]

m6 <- lm(e2_pres_confidence ~ e2_gender * female, data = qualtrics)
temp <- coeftest(m6, vcov = vcovHC(m6, type = "HC1"))
temp
m6$se <- temp[, "Std. Error"]


stargazer(m1, m2, m3, m4, m5, m6,
  dep.var.labels = c(
    "Nom.\\ Confidence",
    "Nom.\\ Confidence",
    "Nom.\\ Confidence",
    "Pres.\\ Confidence",
    "Pres.\\ Confidence",
    "Pres.\\ Confidence"
  ),
  covariate.labels = c(
    "condition: female pronouns",
    "Independent", "Democrat", "female pronouns x Independent",
    "female pronouns * Democrat", "female respondent",
    "female pronouns * female respondent", "(Intercept)"
  ),
  model.names = NULL,
  title = "\\textbf{Effects of Gender Treatment in Nominee Experiment.} The table below shows
          average treatment effects of using female pronouns relative to male in the nominee experiment pooled
          across all doses. Models (a) and (d) show average effects in the entire sample.
          Models (b)-(c) and (e)-(f) condition on respondent party and gender, respectively.",
  se = list(m1$se, m2$se, m3$se, m4$se, m5$se, m6$se),
  label = "table:nominee",
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)



### Table B14  -----------------------------------------
### condition on agency


m1 <- lm(e2_nominee_confidence ~ e2_gender * e2_agency, data = qualtrics)
temp <- coeftest(m1, vcov = vcovHC(m1, type = "HC1"))
temp
m1$se <- temp[, "Std. Error"]

m2 <- lm(e2_pres_confidence ~ e2_gender * e2_agency, data = qualtrics)
temp <- coeftest(m2, vcov = vcovHC(m2, type = "HC1"))
temp
m2$se <- temp[, "Std. Error"]

stargazer(m1, m2,
  dep.var.labels = c(
    "Nom.\\ Confidence",
    "Pres.\\ Confidence"
  ),
  covariate.labels = c(
    "condition: female pronouns",
    "Agency: Education", "Agency: HHS", "Agency: Treasury",
    "female pronouns * Education", "female pronouns * HHS",
    "female pronouns * Treasury", "(Intercept)"
  ),
  model.names = NULL,
  title = "\\textbf{Effects of Gender Treatment in Nominee Experiment by Agency.} The table below shows
          average treatment effects of using female pronouns relative to male in the nominee experiment conditioning on which agency respondents were assigned to read about.",
  se = list(m1$se, m2$se),
  label = "table:nominee_agency",
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)


### Table B15  -----------------------------------------

### Test for Order Effects in the Nominee Experiment


qualtrics$order_exp2 <- factor(qualtrics$order_exp, levels = c("job", "dose", "charts"))

m1 <- lm(e2_nominee_confidence ~ e2_gender * order_exp2, data = qualtrics)
temp1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC1"))
m1$se <- temp1[, "Std. Error"]
temp1

m2 <- lm(e2_pres_confidence ~ e2_gender * order_exp2, data = qualtrics)
temp2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC1"))
m2$se <- temp2[, "Std. Error"]
temp2

stargazer(temp1, temp2,
  column.labels = c("Nom.\\ Confidence", "Pres.\\ Confidence"),
  title = "\\textbf{Effect of Experiment Order on Treatment Effects in Nominee Experiment.} Omitted category is condition presenting nominee experiment first in the survey instrument.",
  se = list(m1$se, m2$se),
  covariate.labels = c(
    "condition: feminine pronouns",
    "framing experiment first",
    "visualization experiment first",
    "feminine pronouns*framing first",
    "feminine pronouns* visualization first",
    "(Intercept)"
  ),
  label = "app:table_exp_order_nominee",
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)






# CHART ANALYSIS  -----------------------------------------



### Table B16  -----------------------------------------


m1 <- lm(e3_serve_US ~ factor(e3_chart), data = qualtrics)
temp <- coeftest(m1, vcov = vcovHC(m1, type = "HC1"))
temp
m1$se <- temp[, "Std. Error"]

m2 <- lm(e3_serve_US ~ factor(e3_chart) * party2, data = qualtrics)
temp <- coeftest(m2, vcov = vcovHC(m2, type = "HC1"))
temp
m2$se <- temp[, "Std. Error"]

m3 <- lm(e3_serve_US ~ factor(e3_chart) * female, data = qualtrics)
temp <- coeftest(m3, vcov = vcovHC(m3, type = "HC1"))
temp
m3$se <- temp[, "Std. Error"]

m4 <- lm(e3_women_rep ~ factor(e3_chart), data = qualtrics)
temp <- coeftest(m4, vcov = vcovHC(m4, type = "HC1"))
temp
m4$se <- temp[, "Std. Error"]

m5 <- lm(e3_women_rep ~ factor(e3_chart) * party2, data = qualtrics)
temp <- coeftest(m5, vcov = vcovHC(m5, type = "HC1"))
temp
m5$se <- temp[, "Std. Error"]

m6 <- lm(e3_women_rep ~ factor(e3_chart) * female, data = qualtrics)
temp <- coeftest(m6, vcov = vcovHC(m6, type = "HC1"))
temp
m6$se <- temp[, "Std. Error"]




stargazer(m1, m2, m3, m4, m5, m6,
  dep.var.labels = c(
    "Government Effective",
    "Enough Women in Government"
  ),
  title = "\\textbf{Effects of Emphasizing Levels over Trends in Government Statistics.} The table below shows
          average treatment effects of presenting data on gender composition in the executive branch in ways that emphasize levels over trends.
          Models (1) and (4) show average effects in the entire sample.
          Models (2)-(3) and (5)-(6) condition on respondent party and gender, respectively.",
  se = list(m1$se, m2$se, m3$se, m4$se, m5$se, m6$se),
  label = "table:charts",
  covariate.labels = c(
    "condition: 2",
    "condition: 3",
    "condition: 4",
    "Independent", "Democrat",
    "Condition 2 * Independent",
    "Condition 3 * Independent",
    "Condition 4 * Independent",
    "Condition 2 * Democrat",
    "Condition 3 * Democrat",
    "Condition 4 * Democrat",
    "female",
    "Condition 2 * female",
    "Condition 3 * female",
    "Condition 4 * female",
    "(Intercept)"
  ),
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)
